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CHAPTER I. INTRODUCTION 


The upper atmosphere (the atmosphere above the troposphere) has 
long been a data-void and therefore, mysterious region. Contrary to 
the lower atmosphere (the troposphere), where easily available data 
contributed to the development of the equations, theories and models 
that form the basis for modem weather forecasting, knowledge of the 
properties of the upper atmosphere were based, until recently, on 
widely scattered wind and temperature measurements derived from meteor 
trails and from balloons penetrating the stratosphere. Now, measurements 
taken regularly by rockets carrying meteorological instruments have 
begun to supply the data which will help solve the mysteries of the 
upper atmosphere. 

From these measurements, in a span of time little more than one 
decade, scientists have compiled a comprehensive climatology of the 
stratosphere and mesosphere. Individual soundings, however, sometimes 
showed large deviations from expected values. Researchers have sought 
to explain these deviations as products of interactions between regions 
of the atmosphere, particularly as consequences of the transfers of 
energy between the lower atmosphere (the troposphere) and the upper 
atmosphere. This thesis examines one such possible interaction 
between the lower and upper atmospheres; namely, the effect of trop- 
ospheric convection in the form of thunderstorms on temperatures 
observed at mesospheric heights. 

From August 15 through September 10, 1976, the National Aeronautics 
and Space Administration (NASA) launched one or two meteorological 
rockets each day from Wallops Island, Va. The rockets reached 
altitudes ranging from 70-90 kilometers and instruments released from 



these rockets measured vertical profiles of temperature, wind, density 
and pressure. Examination of the temperature profiles revealed that 
on three of the days, August 28, September 2 and September 10, a 
portion of the mesosphere (about 60-70 kilometers) exhibited tempera- 
tures 10-15 degrees K wanner than on preceding or following days. 

(See Figure 1 for plots of temperature versus height for two of the 
days, August 28 and September 2). Synoptic records for the three days 
indicate that only on these three days in the period from August 15 
through September 10, intense convective activity in the form of 
thunderstorms and showers occurred at Wallops Island. The storms 
began 4-12 hours before the rocket was launched. Clearly, these 
occurrences make it tempting and also reasonable to try to link the 
thunderstorm activity with the mesospheric heating in some way. 

In this study, atmospheric internal gravity waves are proposed as 
the coupling mechanism. Gravity waves arise when a parcel is displaced 
in a statically stable atmosphere (one in which the environmental 
lapse rate is less than the adiabatic lapse rate) . The restoring 
force produces an oscillation of the parcel with a characteristic 
frequency called the Brunt-Vaisala frequency. The oscillation produces 
a broad spectrum of waves called gravity waves. Among other wave- 
like phenomena in the atmosphere, they fit between the high-frequency, 
short -wavelength acoustic waves and the low-frequency, long-wave- 
length planetary and Rossby waves. Only gravity waves with frequencies 
less than the Brunt-Vaisala frequency are internal gravity waves, i.e., 
propagate energy vertically and only these will be considered. 

The excitation of gravity waves by updrafts penetrating a stable 
layer (common in thunderstorms) , called penetrative convection, has 
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been studied by many researchers (Townsend, 1965, 1966, 1968; Stull, 
1976). Pierce and Coroniti (1966) considered thunderstorms specifi- 
cally as a means of generating gravity waves. Therefore, the fact 
that thunderstorms can be considered a source for vertically propa- 
gating gravity waves is one reason for choosing them as a means of 
coupling the upper and lower atmospheres. 

Another reason is that gravity waves have been observed at upper 
atmospheric heights. Many researchers (Hines, 1960, 1963, 1965; 
Gossard, 1962; Midgely and Liemohn, 1966; Hooke, 1968; Harris et al., 
1969, to name only a few) have documented and developed theories 
accounting for their presence. Also, the temperature profiles for the 
days on which the heating occurred show a superimposed wave-like 
structure with vertical wavelengths ranging from 5 to 10 kilometers 
which is evidence of vertically propagating wave phenomena the size 
of gravity waves. 

A third reason for suspecting that interval gravity waves excited 

by thunderstorms would have an effect on the mesosphere derives from 

the exponential decrease of atmospheric density p with height and the 

principle of the conservation of energy. In order for the wave's 
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kinetic energy per unit volume, pu' , to remain constant as p decreases 
exponentially, the perturbation velocity associated with the wave, u’, 
must increase exponentially. Therefore, waves originally generated 
in the troposphere will have grown considerably by the time they 
reach mesospheric heights and they will have sufficient magnitude to 
affect motions in the upper atmosphere. 

As the gravity waves propagate upwards through the atmosphere, 
kinetic energy will be associated with the velocity perturbations 
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produced by them. The dissipation of this energy can heat the atmo- 
sphere. It is this heating that will be calculated and compared to 
the observed heating in Figure 1. 

Thus, the purpose of this thesis is to account for the observed 
mesospheric heating through the mechanism of vertically propagating 
gravity waves generated by convective activity. 

In Chapter II, the tropospheric sounding for one of the days 
of the observed heating, August 28, 1976, is analyzed to estimate the 
vertical velocity of a thunderstorm updraft on that day. The data 
for August 28 is used throughout this thesis to explain the procedures. 
The same procedures are used on the data of September 2, 1976, and only 
the results are reported. The data of September 10, 1976 , was dis- 
carded since this was the last day of data collection by NASA and it 
is not known for certain that the temperatures returned to normal the 
next day. 

In Chapter III, some basic theory concerning gravity wave propa- 
gation is reviewed, including ray- tracing. Ray- tracing theory assumes 
waves travel in packets or groups that can be identified by a dominant 
wavenumber and frequency. These wave groups travel at local group 
velocities along paths called rays. Whether the ray path of a certain 
wave group reaches the mesosphere depends on its dominant wavenumber 
and frequency as well as the temperature and wind structure of the 
atmosphere. These factors determine if a wave group's energy is 
reflected, absorbed or allowed to propagate. An estimate of the 
spectrum of gravity waves that reach the level of heating is found by 
applying ray— tracing to many combinations of frequencies and wave- 
numbers. 
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In Chapter IV, the results of Chapters II and III along with 
some ideas of Stull (1976) are used to calculate a lower boundary 
condition for the model developed in Chapter V. This model uses 
equations developed by Bretherton (1966) which are solved by a method 
described by Chapman and Lindzen (1970) . 

In Chapter VI, the heating rates are calculated and presented and 
finally, Chapter VII presents conclusions and recommendations for 
future work. 
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CHAPTER II. AN ANALYSIS OF THE SOUNDING OF AUGUST 28, 1976 


A study of internal gravity waves excited by intense tropospheric 
convection should begin with a study of the thunderstorms themselves. 

On August 28, 1976, showers and thunderstorms passed through Wallops 
Island, Va. The first precipitation was recorded at 0700 GMT (3 AM 
local time) but the thunderstorm was strongest from 1100 to 1300 GMT 
(7 to 9 AM local time) during which 1.1 inches of rain fell. No more 
thunder was heard after 1300 GMT (9 AM local time) although the pre- 
cipitation continued to be heavy at times and three quarters of an 
inch of rain fell from then until the showers tapered off around 1700 
GMT (1 PM local time) . The thunderstorm that was reported at Wallops 
Island, Va. was part of a group of thunderstorms and shower activity 
kicked off by a short wave trough aloft passing through the area. They 
formed primarily over the Delmarva penninsula. 

Figure 2 is a vertical plot of the temperature and dewpoint of the 
air above Wallops Island at 1200 GMT (8 AM local time) on August 28, 
1976. The temperatures have been plotted on a skew T-log p diagram. 

The sounding is seen to be quite unstable with a deep layer of moist 
air. This sounding will be used to determine an estimate of an air 
parcel's vertical velocity as it penetrates a stable region and the 
length of time it spends in the stable region before its energy is 
exhausted. 

The level at which moist air becomes saturated by adiabatic ascent, 
the lifting condensation level (LCL) , is found to be 900 millibars on 
this sounding. Above this level, as the saturated parcel travels up 
through the atmosphere, it will follow the saturated adiabat (process 
curve in Figure 2) on the aerological diagram. The saturated adiabat 
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shows how the parcel's temperature changes as it ascends. A saturated 
parcel's temperature will not decrease as fast as the temperature of an 
unsaturated parcel since the process of condensing water vapor into 
a liquid releases heat. The saturated adiabats are constructed on 
the chart with the assumption that all the condensed water vapor falls 
out of the system as soon as it is produced. Also assumed in this 
analysis of the sounding are the restrictions of the parcel method which 
is a method for testing the instability of a parcel of air. The 
assumptions of the parcel method are: 

1. the parcel's movements do not disturb the environment, 

2. the parcel only moves adiabatically , 

3. the pressure of the parcel equals the pressure of the 
environment at every level and 

4. the parcel does not mix with the environment. 

The last assumption is easily recognized as a poor one since only 
one look at the sides and top of a growing cumulus cloud provides 
evidence of turbulent mixing. This mixing brings outside air into the 
parcel in a process known as entrainment. In an attempt to reduce this 
source of error, the graphical correction technique for entrainment 
described by Iribame and Godson (1973) is used. The rate of entrain- 
ment is difficult to determine. It depends on the size of the cloud, 
its stage of development and the velocities of its ascending air 
currents. In this case, it is assumed that at each step in the graphi- 
cal technique, the parcel is diluted by ten percent through mixing 
with external air. 

The sample air parcel begins its journey at the LCL. Above this 
level, the parcel is warmer than the environment, work is done on the 


9 


parcel by buoyancy forces and it accelerates. Eventually, it reaches 
a level at which its temperature equals the temperature of the envi- 
ronment called the equilibrium level (EL) . In this case, the EL is 
9.5 kilometers (285 mb). 

On a skew T-log p diagram, the area between two curves defining 
the state of the environment and the parcel and the initial and final 
pressure levels is proportional to the amount of work done on or by 
the parcel. In this case, between the LCL and the EL, the parcel is 
warmer than the environment and rises, so the environment is doing 
work on the parcel. The area on this diagram is a function of the 
pressure at the EL (285 mb) and the LCL (900 mb) and the mean tempera- 
ture throughout the layer for the parcel (266°K) and the environment 
(265°K). The proportionality constant relating the area and the work 
per unit mass done on the parcel is the universal gas constant 
R (287 Jkg -3 ^" 1 ). The work per unit mass is calculated to be 330 
Jkg It is assumed that all the energy gained by the parcel is 
transformed into kinetic energy and futhermore, that all the kinetic 
energy is associated with the parcel's vertical velocity at the EL. 
With these assumptions, a vertical velocity of 25.7 ms ^ is calculated 
for that parcel. 

Above the equilibrium level, the parcel becomes cooler than the 
environment, the buoyancy forces become negative and the parcel 
decelerates. To calculate the height to which the parcel penetrates 
before it exhausts its energy supply, use is made of an equation 
relating acceleration to the buoyancy of the parcel. 


dw 

dt 



T 

e 


V 


( 2 . 1 ) 
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w is the vertical velocity, g is the acceleration due to gravity, 
is the temperature of the parcel and T e is the temperature of the 
environment. By making use of the fact that dw/dt is equivalent to 
wdw/dz and integrating both sides, an expression for w as a function 
of height z is obtained. 


w^(z) - v^(EL) 


" g 


T - T 
£ 


dz 


’EL 


( 2 . 2 ) 


Using the sounding and successive applications of Sampson's Rule, 
a profile of w with height is calculated and it is discovered that 
this parcel stops rising (w becomes zero) at 12.5 km. With this w 
profile it is possible to integrate the expression dt = dz/w(z) by 
again using Simpson's Rule. The time the parcel is above the equilib- 
rium level is calculated in this manner to be six minutes. 

Both the parcel' s vertical velocity and the length of time it 
spends above the equilibrium level will be used later to compute the 
size of the impulse provided by the convection that excites vertically 
propagating internal gravity waves that in turn heat the mesosphere. 
First, it is necessary to derive some basic relations used to describe 
the propagation of internal gravity waves. 
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CHAPTER III. GRAVITY WAVES AND RAY- TRACING THEORY 

In this chapter, simple equations describing some characteristics 
of internal gravity waves are derived. They are then used to deter- 
mine the paths followed by vertically propagating wave groups. This 
allows the spectrum of gravity waves that reach the level of heating 
to be determined. 

A dispersion equation relating the frequency of a wave to its 
wavenumber can be derived for internal gravity waves. However, some 
simplifying assumptions must be made. To begin with, the atmosphere 
will be considered incompressible. An incompressible atmosphere filters 
out the high-frequency acoustic wave solutions from the wave equation 
and leaves the lower frequency gravity wave solutions. Although an 
incompressible medium may seem unrealistic, equations with this 
simplification have been found to adequately describe gravity wave 
motions observed in the real atmosphere (Gossard and Hooke, 1975). 

This atmosphere will also be isothermal, isentropic, hydrostatic 
and non-rotating. Because of the last condition, the Coriolis force, 
an apparent force appearing in the equations of motion because of the 
earth’s rotation, may be ignored. The Coriolis force is a weak 
deflecting force that is effective only over a long period of time. 

Only periods on the order of less than one hour will be considered here 
so it is reasonable to neglect this force. With these simplifications 
and assuming two-dimensional propagation in the x-z plane, the set of 
equations is: 


du _ _1 3£ 

dt p 9x 


(3.1) 
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dw 

dt 


e . I 3E 

® p 3z 


(3.2) 


dp 

dt 


0 


(3.3) 


3u _3w 
3x 3z 


0 


(3.4) 


For an explanation of all symbols, see the List of Symbols on page vii. 

Equation (3.1) is the equation for the horizontal component of 
motion and Equation (3.2) is the vertical equation of motion. 

Equation (3.3) is an equation representing the incompressibility of 
the fluid and Equation (3.4) is the incompressible continuity equation. 

Equations (3.1) through (3.4) are now linearized by employing the 
well-known perturbation method. Using this method, the dependent 
variables are expressed as the sum of a mean component and a small 
perturbation from that mean. In this case: 

p = p(z) + p' (x,z,t) 
p = p(z) + p’ (x,z,t) 
u = u + u 1 (x, z, t) 
w = w + w* (x, z,t) 

The bar denotes a mean quantity and the prime a perturbation, w is 
assumed to be zero and u, p and p are assumed constant with time t and 
horizontal distance x. u is also constant with height z but p and p 
are allowed to change with height as in an isothermal atmosphere. An 
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important assumption of the perturbation method is that the product 
of two perturbation quantities is assumed to be negligible and may 
be dropped from the equation. Applying the perturbation method to 
Equations (3.1) - (3.4) and letting u = U, a new set of linear equa- 
tions in the perturbation variables is obtained. 



+ U 



. iEl 

9x 


(3.5) 


— r 3w' , TT 3w ’ t 

pt 3T + u &T 1 ■ - gp - 


3z 


(3.6) 


iPl 

3t 



+ w' 



(3.7) 


3u’ 8w T 
3x 3z 


= 0 


(3.8) 


By dividing Equation (3.7) by p and multiplying by g, the third 
term becomes 


i £ i£ = 
V 9 


N 2 w’ 


(3.9) 


N is the symbol for the square of the Brunt-Vaisala frequency. In 

2 — — 

an incompressible atmosphere, N = - (g/p)8p/3z. As stated in the 
introduction, the Brunt-Vaisala frequency is the frequency at which a 
parcel will oscillate if displaced adiabatically an infinitesimal 
distance Az. Thus Equation (3.7) becomes 
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(3.10) 


g 3p f 

•p 3t 




N 2 w’ 


= 0 


Since the motions to be investigated are wave-like, the perturba- 
tion variables will be expressed in the forms 


u' = Uexp[i(cot - kx - mz) ] (3.11) 

w' = Wexp[i(o3t - kx - mz) ] (3.12) 

= Rexp[i(wt - kx - mz) ] (3.13) 

P 

= Pexp[i(wt - kx - mz)] (3.14) 

P 


where only the real parts of the right sides are of interest, k is 
the wavenumber in the x direction and m is the wavenumber in the z 
direction. This form assumes sinusoidal variations in time and in the 
horizontal and vertical directions. Substituting (3.11) - (3.14) into 
equations (3.5), (3.6), (3.10) and (3.8) and noting that the term 
exp[i(wt - kx - mz) ] will be common to all terms and drop out, the 
following are obtained : 


(to - kU)U - kP = 0 (3.15) 

i(w - kU)W + gR - imP = 0 (3.16) 
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(3.17) 


i (w - kU)gR - N 2 W = 0 
kU + mW = 0 . (3.18) 


In matrix form, these equations look like this: 


u-kU 

0 

-k 

0 


U 

0 

i(a3-kU) 

-im 

g 


W 

k m 

0 

0 


P 

2 

0 -N 

0 

i(to-kU) g 


R 


(3.19) 


A trivial solution to this set of equations would beU=W=P=R=0. 
For a non- trivial solution to exist, the determinant formed by the 
coefficients must equal zero. When the equation formed by calculating 
the determinant and setting it equal to zero is solved for the 
frequency a), the result is the dispersion relation: 


a> = kU + N < 


2 2 
k + m 


1/2 


(3.20) 


When two or more waves with slightly differing frequencies and 
wavenumbers propagate at the same time through a medium, they interact 
to form an interference pattern in their amplitudes. These patterns 
are referred to as wave groups or modulation envelopes. The wave groups 
propagate at the group velocity, a speed different from the speed of 
an individual wave. The group velocity can be thought of as the speed 
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at which energy is transmitted by the waves. It is equal to the 
frequency differentiated with respect to the wavenumber. In this case, 
using Equation (3.20), the group velocity in the x direction is: 


3u) 

i ae 

g 3k 


U + N 


m 


2 . 2.3/2 

(k + m ) 


(3.21) 


In the z direction, the group velocity is : 


3oj 

w = + N 

g dm — 


km 


,,2 , 2,3/2 

(k + m ) 


(3.22) 


Bretherton (1966) has shown how equations (3.20), (3.21) and 
(3.22) can be used to calculate the path that a certain wave group 
with a dominant frequency and wavenumber will follow as it moves 
through an atmosphere in which both U and N change with height. For 
these raypaths to be valid, the waves and the atmosphere must obey 
the WKB approx ima tion. That is, the vertical wavelength should be 
small compared to the height over which the horizontal wind U and the 
Brunt-Vaisala frequency N change significantly. In this atmosphere, 

u and w will change with height also. 

8 8 

Bretherton produces the important result that following a wave 
group at its local group velocity, the dominant frequency and horizontal 
wavenumber by which the wave group is identified will remain constant. 

By calculating the vertical wavenumber m from Equation (3.20) at even 
intervals (every kilometer in this case) using local values of U and 
N and a constant w and k, it is possible to determine the vertical 
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profile of m. Equations (3.21) and (3.22) can then be solved for the 
group velocities at each level using local values of m, D and N. The 
time it takes the wave group to reach each new level is computed from 
the local group velocity and the distance between levels (1 kilometer) . 
The distance downwind that the wave group travels can be computed using 
the local horizontal group velocity. In this way, the information 
needed to plot the paths the wave groups travel, their raypaths, is 
obtained. The raypaths are calculated up to 85 km. Figure 3 shows 
the raypaths of several wave groups. 

All the waves propagating upwards through the atmosphere do not 
necessarily reach its uppermost regions. Many are affected by re- 
flection and critical levels. 

If the Doppler-shifted frequency, w - kU, equals the local value 
of the Brunt-Vaisala frequency N, the raypath becomes vertical and 
the group is reflected down towards the ground. An example of this is 
shown in Figure 3. 

A wave group’s critical level is the level at which its Doppler- 
shifted frequency is equal to zero. In the vicinity of a critical level, 
the vertical group velocity becomes proportional to the square of the 
distance from it and the wave takes an infinite time to reach it 
(Bretherton, 1966). Thus, the energy of the wave can't be transmitted 
through this level and is absorbed there. 

Critical levels and reflection are important to consider since 
their existence or non-existence determines which gravity waves can 
reach the mesosphere. This means that the background wind measurements 
are extremely important since they play a role in calculating the 
Doppler-shifted frequency. Whether a critical level exists or not 


18 




Figure 3. Raypaths of gravity wave groups initiated at 15 km 
on August 28, 1976. 
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can depend on a difference in the wind speed of only a few meters 
per second. Therefore, absolute statements cannot be made about the 
spectrum of gravity waves reaching a certain level unless one is 
certain of the accuracy of his wind measurements. However, general 
conclusions may be reached about the size of the spectrum if the paths 
of many waves with different frequencies and wavenumbers are studied. 
This has been done and Figure 4 is a plot of the fate of gravity waves 
of different frequencies and wavenumbers (plotted as periods and 
horizontal wavelengths) in the atmosphere of August 28, 1976. This 
plot will be used in the next chapter to determine the range of 
frequencies and wavenumbers that can penetrate this atmosphere to the 
level of heating. 
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Figure 4. Fate of gravity wave groups identified by dominant wavelength 
and periods in atmosphere of August 28, 1976. 
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CHAPTER IV. CALCULATION OF THE CONVECTIVE ENERGY INPUT 


In Chapter II, it was shown that a parcel rising in unstable air 
reaches an equilibrium level and then penetrates a distance into a 
stable layer, before sinking back down. This type of behavior has been 
studied in the context of thermals overshooting the atmospheric 
boundary layer into an overlying temperature inversion and has been 
labeled penetrative convection. Stull (1976) has studied this type of 
penetrative convection and many of his ideas can be applied to describe 
the penetrating elements found in a thunderstorm. 

Stull has developed a useful representation of an idealized 
convective element. See Figure 5. A parcel penetrates a distance d 
into a stable region and then sinks back to the equilibrium level. 

Figure 5a represents the parcel’s height above the equilibrium level 
as a sinusoidal function of time. P/2 is the time the parcel is above 
the equilibrium level, calculated in Chapter II to be six minutes. 

Figure 5b is the idealized element at the time of maximum penetration 
as a function of distance x. The figure shows that some of the stable 
air is drawn down below the equilibrium level to conserve mass. D is 
the diameter of the element at the equilibrium level. 

In the real atmosphere, thunderstorm updrafts acting as penetra- 
tive elements excite internal gravity waves with a broad range of 
frequencies and wavenumbers. The irregularities of the individual 
elements make the precise determination of the frequencies and wave- 
numbers difficult. Stull solves this problem by making a harmonic 
analysis of his idealized disturbance. He comes up with a smooth, almost 
Gaussian frequency spectrum and an equally simple shape for the wave- 
number spectrum. Both are shown in Figure 6. 
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Note from the frequency spectrum that the bulk of the energy is 
associated with low-frequency or longer period waves. From the wave- 
number spectrum it is not as easy to see with which wavenumbers most 
of the energy is associated but there is a wavelength at which the 
energy density is a maximum and that may be calculated. 

Stull combines the equations describing the shapes in Figure 6 
to obtain an expression for the spectral energy density J of one 
idealized convective element. 


J = 


’'H 


2 in ?'W 

II 


exp 



(4.1) 


where f = 
is the 
number in 


2u/P and = tt/D. uk is the Doppler-shifted frequency and 

2 2 2 

horizontal wavenumber (k^ = k + £ where k(£) is the wave- 

the x(y) direction) . J is normalized so that 


r 2 tt t 00 


OO 


JdUhdk^d(J> = 


= 1 


4>=0 kg =0 w -0 


(4.2) 


(J) is the horizontal angle in wavenumber space. 

If the normalized energy density function (Equation (4.2)) is 
integrated over a range of wavenumbers and frequencies corresponding 
to the wavelengths and periods of the range of waves calculated to 
reach the level of heating by ray- tracing theory, the result will be a 
rough estimate of the fraction of the total kinetic energy of the 
parcel that is used to excite those gravity waves at the equilibrium level. 
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The mean value theorem is used to integrate the expression 

JdULdk^d<j> = JAw^Ak^Acj) . (A. 3) 

$ kg ^ 

J is a mean value of the energy density function, A<*h is the frequency 
interval and Ak^ is the horizontal wavenumber interval. Since the 
energy density function J is integrated over all horizontal angles in 
wavenumber space, A<}> equals 2 tt. The next step is to calculate 
and Aco^ with the help of Figures 4 and 6. 

The wavenumber at which the energy density function is a maximum 
(the peak of the graph in Figure 6b) is calculated from Equation (4.1) 
to be = tt/D. In terms of horizontal 
wavelength at maximum J is equal to 2D. Assuming the diameter of the 
element at the equilibrium level is one kilometer, the wavelength for 
maximum J is two kilometers. Combining this fact with the information 
gathered from Figure 4 that wavelengths less than twelve kilometers do 
not reach the level of heating and noting that wavelength increases 
toward the origin in Figure 6b, it can be seen that the range of wave- 
lengths contributing the energy to the heating, Akg, will be found to 
the left of the peak and that twelve kilometers is the right edge (or 
short-wavelength side) of the range. 

The chart in Figure 4 shows that the long-wavelength side of the 
wavelength range is not easily defined. All the longer wavelengths com- 
bined with the range of frequencies that can be considered reach the 
level of heating , although Figure 6b shows that little energy is 
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associated with the very long wavelengths. There are two reasons for 
no cut-off. First, the larger the wavelength, the smaller the wave- 
number and therefore the Doppler-shifted frequency (uHk.U) will never 
be zero (with the range of frequencies considered here). Thus, a 
critical level will never be reached. Also, a small wavenumber means 
that even with the strong easterly winds of the stratosphere, the 
Doppler-shifted frequency, ok, will always be less than the Brunt- 
Vaisala frequency and the wave group will not be reflected. 

Since no cut-off on the small- wavenumber side of Figure 6b is 
indicated from ray-tracing theory, the wavenumber interval is estimated 
to range from the clearly defined edge at the wavenumber associated 
with a wavelength of 12 kilometers all the way down to the origin 
(wavenumber equal to zero). Therefore, Ak^ is approximately 0.5 km 
The frequency interval. Aw , can be estimated from Figure 4 by 
looking at the area on the figure associated with the most energy for 
the wavelength range and assuming the frequency range there. This 
range taken at wavelength equal to 12 kilometers is 0.0013 sec \ 

The fractional part of the energy density associated with the 
range of gravity waves that reach the mesosphere is found by calculating 
the right side of Equation (4.3) using the values for AW. and Ak just 
estimated. The fraction is 0.0021. If the total energy associated 
with the rising parcel calculated in Chapter II from the temperature 
sounding of August 28 is multiplied by this fraction, the result is 
the amount of energy associated with the range of gravity waves that 
reach the level of heating. The assumptions that this energy is all 
converted into kinetic energy and all the kinetic energy is associated 
with a vertical velocity are made. Therefore, 
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KE = ~ w 2 (4.4) 

where KE represents the fractional kinetic energy per unit mass 
(0.7Jkg ^) . From Equation (4.4), a vertical velocity of 1.2 ms ^ is 
calculated. 

'i 

The kinetic energy associated with a vertical velocity of 1.2 ms 
represents an estimate of the energy initially contained by those waves 
calculated to reach the mesosphere by ray- tracing theory. By using this 
value as the impulse that generates the internal gravity waves in a 
model, the energy attenuation by reflection and critical level inter- 
action of waves on their way to the mesosphere will be accounted for. 
Therefore the calculated heating will be representative of the actual 
energy involved. 
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CHAPTER V. THE MODEL OF GRAVITY WAVE STRUCTURE 

The equations used to find the detailed gravity wave structure 
are taken from Bretherton (1966) . They are derived assuming a com- 
pressible atmosphere, but the equations describing the propagation of 
wave groups derived previously for an incompressible atmosphere may be 
used in conjunction with them as long as the WKB approximation is 
valid. Also, the Brunt-Vaisala frequency is now expressed as a function 
of the change of potential temperature with height. 

N 2 = f£ 


The linearized equations corresponding to such an atmosphere are 
reproduced from Bretherton (1966) . 
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x or z used as a subscript denotes partial differentiation with respect 

D 3 B 

to that variable. — stands for the operator t— ■ + U - 5 — and y = c /c . 

Dt a t dx p v 

In these equations, p is a local value of the background density 
and is a function of height. It is the exponential variation of p that 
allows the amplitude of the waves to grow with height as a result of 
the conservation of energy. Bretherton rescales the equations to 
account for this effect by naming a reference density p. He defines 
new variables in this way: 

(u* > w’) = (p/p ) -1/2 (u,w) (5.6) 
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The equations used to find the detailed gravity wave structure 
are taken from Bretherton (1966) . They are derived assuming a com- 
pressible atmosphere, but the equations describing the propagation of 
wave groups derived previously for an incompressible atmosphere may be 
used in conjunction with them as long as the WKB approximation is 
valid. Also, the Brunt-Vaisala frequency is now expressed as a function 
of the change of potential temperature with height. 
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The linearized equations corresponding to such an atmosphere are 
reproduced from Bretherton (1966) . 
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x or z used as a subscript denotes partial differentiation with respect 

D 3 3 

to that variable. — stands for the operator w— + U x— and Y = c /c . 

Dt r dt dx 1 p v 

In these equations, p is a local value of the background density 
and is a function of height. It is the exponential variation of p that 
allows the amplitude of the waves to grow with height as a result of 
the conservation of energy. Bretherton rescales the equations to 
account for this effect by naming a reference density p. He defines 
new variables in this way: 
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Equations (5.2) to (5.5) become 


u + U w + — p = 0 
Dt z x 


(5.9) 


2- 5 + g+ I c |_-I)p, 0 


(5.10) 


1 1 D ~ . ~ . ,3 . l.~ „ 

^=M p + D k + < 3^ + H ) ' ,= 0 
c p 


(5.11) 


D ~ 2- 

— CT - N w = 0 


(5.12) 


30 



Here, c is the speed of sound equal to yp/p and H is a measure of the 
scale height such that: 


1 

H 


2p P z + 0 6 z 


(5.13) 


At this point, c is allowed to approach infinity as in an incompress- 
ible atmosphere with the sole effect being the elimination of acoustic 
wave solutions. This leaves the continuity equation (Eq. 5.11) in the 
form 


U x + ( -fe “ H> * = ° • (5 ' 14) 

(5.9), 


(5.15) 

The first term on the right-hand side comes from the approximation 
Bretherton makes that w = - w /H, 

XZ X 

By assuming a wave form in the horizontal direction and time with 
a real wavenumber and frequency, the vertical structure of the wave may 
be found. Therefore, 

w = w(z) exp[i(kx - u)t) ] 


As a result of eliminating u, p, and d among Equations 
(5.10), (5.12) and (5.14), the following equation is found: 
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and the equation becomes: 
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This equation is of the form 

w + f(z)w = 0 (5.17) 

zz 

where f(z) refers to the terms within the brackets and is a function of 
height. 

Although this equation resembles a simple wave equqtion, the 
function of height f(z) prevents it from being solved analytically. 
However, a numerical technique described by Chapman and Lindzen (1970) 
and originally discussed by Bruce, Peaceman, Rachford and Rice (1953) 
may be used. This method uses finite difference approximations for 
the derivatives and incorporates boundary conditions for the top and 
bottom of the model. The bottom boundary condition is a vertical 
velocity of 1.2 ms ^ at a level of ten kilometers. Ten kilometers is 
the approximate height of the equilibrium level for a penetrative 
element on August 28, 1976. A radiation condition is used as the upper 
boundary condition. It assumes that there are no energy sources at 
z - °°. 

The model’s domain extends from 10 km to 80 km and the region is 
divided into 100-meter- thick layers. This fine resolution is necessary 
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for numerical reasons (see Chapman and Lindzen, 1970) . All the data 
needed to solve the equation (winds, temperature, scale height, etc.) 
are calculated or taken directly from the rocket measurements. Linear 
interpolation is used to find values between measurement levels. 

The upper level winds are a smoothed profile from wind measurements over 
several days. See Figure 7. 

The equation is solved with a predetermined frequency and horizontal 
wavenumber. The values for to and k are only those combinations found 
to describe wave groups that reach the level of heating. 

The resulting values of w are multiplied by the factor in 
Equation (5.6) to put back the density effect. Values of w’ , the actual 
vertical velocity perturbations resulting from the propagation of 
internal gravity waves of some frequency u) and wavenumber k, are then 
obtained. A plot of the amplitude of w' (since w’ is complex) is 
shown in Figure 8. It reveals the great effect decreasing density with 
height has on the structure of the waves. 

The horizontal velocity perturbation u’ can now be calculated 
using the continuity equation (Eq. 5.14), since the vertical variation 
of w* and thus w is now known. A wave form is assumed for u in x and 
t, i.e., k and to are made real, so that 

u * u(z) exp[i(kx - tot) 3 


and Equation (5.14) becomes 


i ,3 1. ~ 

u * f ai * h> w 


(5.18) 
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Solving this equation for u and multiplying the result by the factor 
in Equation (5.6), the horizontal perturbation velocities corresponding 
to the internal gravity waves are obtained. A plot of the amplitude 
of u T is shown in Figure 9 in which the density effect on wave 
amplitudes is also apparent. 
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AMPLITUDE OF VERTICAL PERTURBATION VELOCITY (M S"') 


Figure 8. Amplitude of vertical perturbation velocity w' as a 
function of height. 
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Figure 9. Amplitude of horizontal perturbation velocity u* as a 
function of height. 
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CHAPTER VI. CALCULATION OF HEATING RATES 


The gravity wave solutions found in Chapter V were obtained 
without any friction. The rate at which friction would dissipate 
energy in those solutions is now calculated and it is assumed that 
this energy heats the background atmosphere. 

Gravity waves have been shown to have the potential to produce 
large amounts of heating. Hines (1965) has studied this effect at 
heights above the mesopause (80 km) . He suggests that heating by 
gravity waves may be responsible for the rapid rise of temperature 
with height that occurs in the E-region (100-120 km) . Whitaker (1963) 
has suggested that gravity waves continuously propagating upwards 
from the photosphere of the sun are responsible for maintaining the 
high temperature of the solar corona. 

In this case, an estimate of just how much of the kinetic energy 
of these perturbations is destroyed by frictional dissipation is found 
by assuming that the friction term may be represented by means of an 
eddy viscosity coefficient K. K is a function of height. Cunnold, 
Alyea, Phillips and Prinn (1975) developed the profile of K for the 
stratosphere and mesosphere seen in Figure 10 which will be used here. 
In this case, the wind varies only in the vertical and the viscous 
dissipation terms can be represented this way: 
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In order to determine the rate of destruction of kinetic energy 
(KE) by these terms, they are multiplied by u’ and w’ respectively 
and added. 


.,2 
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dt 
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dt dt 3 Z 3z ' 
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W 3z (K 3z ' 


( 6 . 1 ) 


Next the thermodynamic equation may be employed'. 


dq dT . da 

—a. - c ___ + p — - 

dt v dt r dt 


( 6 . 2 ) 


c^, a constant, is the specific heat at constant volume and a is the 
specific volume. The kinetic energy destruction rate is set equal to 
the heat added per unit time, dq/dt, since it is assumed that all the 
kinetic energy due to the velocity perturbations caused by gravity 
waves is dissipated as heat. 

In order to obtain a heating rate from this equation, the assump- 
tion must be made that all the heat energy per unit time, dq/dt, goes 
into heating the local medium as represented by the term c^dT/dt and 
none into expanding the medium as represented by the term pda/dt, since 
it is difficult to estimate the amount of heat used in expansion. 
Therefore, it must be remembered that the heating rates obtained this 
way (in the form degrees K per day) will be the upper limit of the 
actual heating possible. 
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Some results obtained by solving 


dT _ 1_ d^i 

dt “ c dt 
v 


(6.3) 


are shown in Figures 11, 12 and 13. 

All the procedures performed on the data of August 28 are now 
used with the data of September 2. Convective activity occurred on 
that date as a result of a cold front moving through the area. 

When raypaths are calculated for wave groups in the atmosphere 
of September 2, some changes are necessary in the low level wind 
profile. The sounding for that day at 8 AM local time recorded 
westerly winds on the order of 30 m/ sec at the levels where the waves 
would be initiated by the penetrative convection. These strong 
positive winds produce critical levels (w - kU = 0) for a large range 
of frequencies and wavenumbers, thus effectively putting a lid on the 
propagation of most of the wave energy. 

Since a sounding only reflects the value of the wind speed at one 
instant and at one point in space and since penetrative convection was 
going on all around the area before and after the time of the sounding, 
it is not unreasonable to assume that the wind can have different 
values at other times. Therefore, the wind values used in computation 
of the raypaths have been changed to smaller positive values (y 5 ms 
to allow more of the wave groups to reach the mesosphere. 

Some heating rates calculated for September 2 may be seen in 
Figures 14, 15 and 16. 
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Figure 11. Heating rates as a function of height, August 28, 1976. 
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Figure 12. Heating rates as a function of height, August 28, 1976. 
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Figure 13. Heating rates as a function of height, August 28, 1976. 
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Figure 14. Heating rates as a function of height, September 2, 1976. 
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Figure 15. Heating rates as a function of height, September 2, 1976. 
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Figure 16. Heating rates as a function of height, September 2, 1976. 
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CHAPTER VII. CONCLUSIONS 


The purpose of this thesis was to determine whether heating 
produced by gravity waves at an altitude of 60-70 kilometers could 
be responsible for the rapid heating observed in rocket measurements. 

The results seem to indicate that gravity waves are the probable 
mechanism producing that heating, despite all the assumptions and 
approximations made. 

A review of the assumptions is presented here. 

1. The assumptions of the parcel method (listed in Chapter II) 

2. The assumption that all the parcel’s kinetic energy becomes 
associated with the parcel's vertical velocity at the EL 

3. The assumptions associated with ray-tracing theory 

a. the background atmosphere is 

incompressible 

isothermal 

isentropic 

hydrostatic 

non-rotating 

b. waves and atmosphere must obey the WKB approximation, 
i.e., variables are slowly varying over a wavelength 

4. The assumptions associated with the set of equations used to 
find the vertical gravity wave structure in Chapter V 

a. again, the WKB approximation 

b. nonlinear terms are small and may be ignored 

c. a wave form for w and u with real wavenumber and frequency 

d. compressible and non-rotating atmosphere 

e. no friction 

5. The assumption that the amount of frictional dissipation is 
calculated using waves that traveled to mesospheric heights with 
no friction 

6. The assumption that the friction term may be estimated by 
means of an eddy viscosity coefficient 
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7. The assumption that all the energy dissipated by friction 
goes into heating the local medium 
From this list of assumptions, the limitations of this particular 
approach of calculating the heating rates become readily apparent. 

The large amplitudes of u’ and w’ evident particularly above 70 
kilometers in Figures 8 and 9 make it unreasonable to assume linear 
solutions as in this model. Therefore, rates obtained above 70 kilo- 
meters are not considered valid. Even below 70 kilometers the 
amplitudes in some cases are large enough to demand that nonlinear terms 
be accounted for. The large amplitudes and the complications intro- 
duced by the nonlinearities may be reduced if a friction term is 
included in the equations of motion used to determine the detailed 
structure of the gravity waves. This was not done in this study 
because of the great difficulty of solving such a set of equations. 

Most of these assumptions (particularly the last three) will lead 
to overestimation of the heating rates. Therefore, large rates will 
be expected but if some correction factor could be applied to the rates 
accounting for all the overestimations, heating similar to that ob- 
served would very likely be the result. The large rates shown in the 
figures then tend to substantiate the claim that the heating is indeed 
a product of gravity wave interaction with the atmosphere. 

The waves also produce cooling in certain regions. This is easily 
explained. Dissipation by friction is a process in which perturbations 
in a system become smoothed and flattened out. This is done not only 
by depositing energy in regions (heating) but also by removing energy 
from other regions (cooling) . Note that even with cooling there is al- 
ways net heating in the 60-70 kilometer range. 
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Another interesting feature in the plots of the heating rates is 
revealed through the shape of the plots and not the actual values. In 
all the plots, the calculated heating begins to take on significant 
values at the level at which the heating was actually observed. It 
seems as if heating by gravity waves only takes on significance above 
approximately 60 kilometers which may be why the heating is only 
noticeable above that level. 

If gravity waves generated by convection in thunderstorms do 
produce heating of the magnitude observed consistently (other processes 
in the atmosphere may mask the effect at times) , then in regions of the 
earth where thunderstorms and intense convection occur frequently, such 
as the tropics, the mesosphere may have slightly higher temperatures. 

This remains to be determined. 

This study has produced no evidence that gravity waves are not the 
mechanism responsible for the observed heating. Any mechanism capable 
of coupling two regions of the atmosphere is potentially important so 
it is recommended that now a more detailed study should be made involving 
fewer simplifying assumptions and using equations for the gravity wave 
structure that include a friction term. 

Since intense convection occurs somewhere on the earth practically 
every minute of every day, upper atmospheric measurements made at 
other sites should eventually coincide with the passage of a squall 
line. Such occurrences should be watched for. Similar warmings in 
other locations will validate the Wallops Island data and provide more 
cases for study. A special program involving a series of rocket 
launchings before, during and after intense convective events would also 
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help to better determine the existence, extent and intensity of the 
heating. 

In Chapter III, the importance of the background wind in deter- 
mining the paths and range of gravity wave groups was discussed. 

Given that importance, a study of the variability of stratospheric 
and mesospheric winds would be of use. A look at the sensitivity of 
wave propagation to wind and temperature Values would also be of 
interest. 

All of these suggestions will lead to a much-needed understanding 
of how events occurring in distant and seemingly distinct regions of 
the atmosphere are related. It is hoped that this understanding will 
serve as another step towards understanding the behavior of the atmo- 
sphere as a whole. 
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